#include <algorithm>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <iterator>
#include <list>
#include <map>
#include <set>
#include <sstream>
#include <string>
#include <queue>
#include <vector>
#include <cmath>
#include <boost/program_options.hpp>
#include "MetastableSet.h"
#include "StringSplit.h"
#include "BlockAverage.h"

#define MaxLineLength 2048

namespace po = boost::program_options;
using namespace std;

inline int angleIdx (const double & myangle,
		     const double & bin)
{
  return int((myangle + 180) / bin);
}


int main(int argc, char * argv[])
{
  std::string ifile, ofile;
  double center, width, bin;
  unsigned nBlock;

  po::options_description desc ("Allow options");
  desc.add_options()
      ("help,h", "print this message")
      ("meta-center,c",  po::value<double > (&center)->default_value (180),   "center of the metastable set")
      ("meta-width,w",  po::value<double > (&width)->default_value (15),   "width of the metastable set")
      ("bin,r",  po::value<double > (&bin)->default_value (5.),   "the bin size")
      ("num-block,n",  po::value<unsigned > (&nBlock)->default_value (10),   "number of block")
      ("input,f", po::value<std::string > (&ifile)->default_value ("angaver.xvg"), "input angle file")
      ("output,o", po::value<std::string > (&ofile)->default_value ("fmpt.out"), "the output of first mean passage time");
  
      
  po::variables_map vm;
  po::store(po::parse_command_line(argc, argv, desc), vm);
  po::notify (vm);
  if (vm.count("help")){
    std::cout << desc<< "\n";
    return 0;
  }

  unsigned nbin = (360. + 0.01 * bin) / bin;
  bin = 360. / double(nbin);
  cout << "# bin size: " << bin << endl;
  cout << "# nbin: " << nbin << endl;
  // record
  vector<vector<double > > record (nbin);

  double metal, metau;
  metal = center - width;
  if (metal < -180) metal += 360.;
  int tmpn = (metal + 180 + 0.5 * bin) / bin;
  metal = -180 + tmpn * bin;
  metau = center + width;
  if (metau > 180) metau -= 360.;
  tmpn = (metau + 180 + 0.5 * bin) / bin;
  metau = -180 + tmpn * bin;
  cout << "# metastable set: [ " << metal << " , " << metau << " ]" << endl;
  // metastable set
  MetastableSet ms (metal, metau);

  vector<double > times;
  vector<double > anglev;
  ifstream fpname (ifile.c_str());
  if (!fpname){
    std::cerr << "cannot open file " << ifile << std::endl;
    return 1;
  }
  char nameline [MaxLineLength];
  while (fpname.getline(nameline, MaxLineLength)){
    if (nameline[0] == '#' || nameline[0] == '@') continue;
    vector<string > words;
    StringOperation::split (string(nameline), words);
    if (words.size() < 2) {
      cerr << "wrong file format of " << ifile << endl;
      exit (1);
    }
    times .push_back (atof(words[0].c_str()));
    anglev.push_back (atof(words[1].c_str()));
  }
  
  unsigned nextstab;
  nextstab = 0;
  for (unsigned ii = 0; ii < times.size(); ++ii){
    if (ms.inSet(anglev[ii])) {
      // record[angleIdx(anglev[ii], bin)].push_back (0.0);
    }
    else {
      if (nextstab > ii){
	record[angleIdx(anglev[ii], bin)].push_back (-(times[ii] - times[nextstab]));
      }
      else {
	nextstab = ii;
	do {
	  nextstab ++;
	  if (nextstab >= anglev.size()){
	    goto line100;
	  }
	} while (! ms.inSet(anglev[nextstab]));
	record[angleIdx(anglev[ii], bin)].push_back (-(times[ii] - times[nextstab]));
      }
    }
  }
  line100:

  BlockAverage ba;  
  FILE * fout = fopen (ofile.c_str(), "w");
  if (fout == NULL){
    cout << "cannot open file " << ofile<< endl;
    exit (1);
  }
  for (unsigned ii = 0; ii < record.size(); ++ii){
    double myangle = -180 + bin * (ii + 0.5);
    if (record[ii].size() == 0){
      fprintf (fout, "%f\n", myangle);
    }
    else if (record[ii].size() < nBlock * 5){
      ba.processData (record[ii], record[ii].size());
      fprintf (fout, "%f %e %e\n", myangle, ba.getAvg(), ba.getAvgError());
    }
    else {
      ba.processData (record[ii], nBlock);
      fprintf (fout, "%f %e %e\n", myangle, ba.getAvg(), ba.getAvgError());
    }
  }  
  fclose (fout);
  
  return 0;
}
